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5 TO WHOM IT MAY CONCERN, THE FOLLOWING IS 

A SPECIFICATION OF THE AFORESAID INVENTION 



ENDPLATE DETECTION IN DIGITAL RADIOGRAPHY BY 
DYNAMIC PROGRAMMING USING BOTH LOCAL AND 
GLOBAL CONSTRAINTS 

CROSS REFERENCE TO RELATED APPLICATION 

The present application claims priority to United States Provisional Patent Application 
Serial No. 60/326,358, filed October 1, 2001, which is hereby incorporated by reference. 

BACKGROUND OF THE INVENTION 

1. Field of the Invention 

The present invention relates generally to analysis of X-ray images, and more 
particularly, to the automatic detection of vertebra endplates for measuring deformities of 
pathological spines in digital radiography. 

2. Description of Related Art: 

In a previously disclosed invention, "Detection of vertebra endplates in digital 
radiography", Attorney Docket No. 00P7820US, an evidence reasoning approach is proposed. 
First, local pieces of evidence of different kinds are collected. They are then combined by a rule- 
based deduction to incorporate certain geometric constraints about endplate shape to determine 
the positions of endplates. Because the reasoning about endplate positions in the previously 
disclosed technique is mainly local, it is difficult to apply some global shape constraints. Global 
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constraints are especially useful when there are other organs interfering with the endplate, 
causing part of the endplate image to be blurred, obscured, or underexposed. Errors in local 
reasoning also tend to propagate: if a mistake is made at one place, it will affect the decision 
making at other places. 

In this invention, a global reasoning approach is adopted. Decision-making about 
endplate positions at local regions is deferred until information is gathered about all endplates in 
the image. This is achieved via a dynamic programming approach. Dynamic programming is a 
general method for optimization. The invention formulates the endplate detection problem in a 
form suitable to solve by dynamic programming. Besides, the conventional dynamic 
programming approach is generalized to allow global constraints to be enforced in a progressive 
manner. This is realized by introducing backward tracing during forward computation of 
endplate scores. As a global constraint on vertebra shape, an analytic model about vertebra 
height distribution is proposed. 



SUMMARY OF THE INVENTION 

Disclosed is a method and apparatus for detecting endplates of vertebra, comprising the 
steps of providing a filtered curvature map derived from filtering an intensity curvature map of a 
spine image in a direction relative to a spine axis, computing minimum and maximum curvature 
projections for each point along said spine axis, at a plurality of points on said spine axis, 
computing a score from said curvature projections, said score indicating the likelihood that the 
point is located on an endplate, and subjecting said score calculation to at least one of the 
following constraints (1) that the angle between neighboring endplates not exceed a certain 
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value, and (2) that the variation in height of vertebra of said spine image should satisfy a global 
model. 

In another aspect of the method said step of computing minimum and maximum 
curvature projections further comprises the steps of computing said spine axis from a spine 
image boundary, at each point p along said spinal axis, compute a curve set S p of all fitting 
curves passing through point p from one boundary to the other that fit the parabolic equation y = 
kx 2 + b, deriving said maximum curvature projection as the maximum of the curvature sums 
along said fitting curves, and deriving said minimum curvature as the minimum of the curvature 
sums along said fitting curves. 

In another aspect of the method said score E is computed from the equation 

N V 

E = 2 C SiWaK (e t ) + ]T max {C 5>min (p) } 

where V is the number of vertebrae, the first term is the sum of the values of said maximum 
curvature projections at endplate positions, and the second term is the sum of the values of said 
minimum curvature projections at intervertebral disc positions. 

In another aspect of the method said score is computed by dynamic programming. 

In another aspect of the method said dynamic programming comprises the steps of 
obtaining the range of first and last endplates, using the range of said first endplate and said 
constraints to determine the range of the next endplate if still within range of said last endplate, 
backward tracing through the state space in the dynamic programming optimization of the score 
to derive a vertebra height profile, fitting a global endplate height model to an augmented height 
profile up to a threshold value of the fitting residual, repeating said step of using said constraints 
to determine the range of endplates on subsequent endplates until the range of the last endplate 
is reached. 



BRIEF DESCRIPTION OF THE DRAWINGS 



Figure 1 is a flowchart diagram of an embodiment of the invention. 
Figure 2(a) is a typical spinal image. 

Figure 2(b) is the same image as 2(a) with superimposed anatomical markings 
Figure 3(a) is a filtered curvature map of Figure 2(a) 
Figure 3(b) shows a curvature projection. 

Figure 4 are graphs of maximum and minimum curve projections. 
Figure 5 is a diagram showing endplate detection as a path finding problem. 
Figure 6 is a flowchart of a method of dynamic programming. 
Figure 7 shows a spinal image before and after endplate detection. 



DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 

It is assumed in this invention that the spine boundary has been detected by a previously 
disclosed technique ("Automatic detection of spine axis and spin boundary in digital 
radiography", Attorney Docket No. 00P7817US), the disclosures of which are incorporated by 
reference herein in their entirety. Based on the spine boundary information, the curvature 
extreme map is extracted and filtered, in the same way as was disclosed in 00P7820US, the 
disclosures of which are incorporated by reference herein in their entirety, but with one 
important improvement: projections are based on a curved endplate model. The curvature map 
projection produces two curves: the maximum curvature projection and the minimum curvature 
projection. These two projections are combined to derive the endplate positions by a dynamic 



programming approach that can consider both local and global constraints on vertebra shape. 
Figure 1 shows the workflow of the detection. 

Referring to Figure 1, control blocks that are identical to those disclosed in 00P7820US 
are labeled with double-digit numerals, while those control blocks new to this invention are 
5 labeled with triple-digit numerals. It should be noted, however, that while blocks 20 and 22 are 
identical to that disclosed in 00P7820US, the method by which these two blocks are executed is 
novel to this invention. At block 10, a spine image is provided, in which endplate orientations 
are desired. This image of a spine may have been taken by X-ray or other imaging technologies, 
such as computerized axial tomography (e.g., CAT scan), sonogram, magnetic resonance (MRI) 

1#0 or other techniques. The image is preferably converted to or taken in digital form. Figure 2 is an 

O 

example of a spine image, which may be provided as input to the endplate detection system of 

m 

j: the present invention. 

u 

m In block 12, a region-of-interest (ROI) is defined, for example, based on boundaries of 

p the spine, which are input from block 100. The detection is preferably confined to the spine 

v, a 

PLJ5 boundary and may be performed by analyzing the image of the spine. 

O In block 110 a curvature map of the spine is extracted, such as by the method disclosed in 

00P7817US. The map is shown in Figure 3(a). 

In blocks 20 and 22 maximum and minimum curvature projections are derived as more 
fully explained below. 

20 In block 130, the dynamic programming disclosed in this invention is executed as will be 

more fully described below. 

In block 30, the best scores provide the endplate positions and orientations. 
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For notational convenience, Figure 2a shows a typical spine image. Figure 2b shows the 
same image labeled with the anatomical landmarks relevant to this invention. 

In the following, two major components of this invention are disclosed in details. They 
are a) the improved curvature projection, and b) the dynamic programming approach for deriving 
endplate positions. 

Curved curvature map projection 

Curvature map projection is performed for the purpose of computing of the local score of 
support for endplates and intervertebral discs. Curvature map projection operates on a filtered 
curvature map. Figure 3(a) shows the filtered curvature map of Figure 2(a), obtained by the 
method disclosed in 00P7820US or other suitable method. 

Figure 3(b) illustrates how a curvature map is projected to obtain curved curvature-map 
projections. 

First, the spine axis is computed from the spine boundary. For each point p on the spine 
axis, draw two lines at a fixed angle to the orthogonal line of the tangential line of the spine axis 
at p . These two lines intersect with the left and right spine boundaries at L x 9 L 2 and 

R x R 2 respectively. The segment L x ,i 2 determines the range of the left end-points of endplates 
on the left spine boundary, while R x R 2 determines the range of the right end-points of endplates 
on the right spine boundary. For each point p L on the curved segment L x L 2 , p R on segment 
R x R 2 , fit a parabolic curve: 



y = kx 2 +b 



to the three points, p L , p 9 and p R . Denote this curve by C (p L9 p 9 p R ). If there are n L points 
on L x L 2 and n R points on R x R 2 , the total number of fitting curves at p is n L • Denote this 
curve set by S p ={c(p i} ) }. 

Referring to Figure 1, box 20, for the maximum curvature projection, the curve set 
S p represents all possible endplates passing through point p . The maximum curvature 

projection at p is obtained as the maximum of the curvature sums along such fitting curves. If 
there were an endplate passing through p , this endplate will be captured by the curve producing 
the maximum curvature projection, because endplates appear as ridges in the curvature map. 
The orientation of the endplate can be directly derived from the fitted curve. 

Referring to box 22, for the minimum curvature projection, the curve set S p represents 

all possible intervertebral discs passing through point p. The minimum curvature projection at 

p is obtained as the minimum of the curvature sums along all these fitting curves. If there were 
an intervertebral disc at p , its position will be captured by the curve producing the minimum 
curvature projection, because an intervertebral disc usually appear as a valley in the curvature 
map. 

The curvature projections are computed for each point on the spine axis. So, two curves 
are produced, one for the maximum curvature projection, denoted by (p) , and one of the 
minimum curvature projection, denoted by (p) . 

Figure 4 (a) and (b) show, respectively, the maximum curvature projection and the 
inverse of the minimum curvature projection of Figure 3 (a) in graphs plotting pixels versus the 
curvature value. The inversion for the minimum curvature projection is intended for just 



converting valleys to peaks. Therefore, peaks on the maximum curvature projection indicate 
possible positions of endplates on the spine axis, whereas peaks on the curve of minimum 
curvature projection mark possible positions of intervertebral disc on the spine axis. 

Dynamic programming to determine endplate positions: 

According to the characteristics of the curvature projection profiles, an endplate should 
correspond to a peak on the maximum curvature projection, and between two endplates - one 
corresponding to the upper endplate of a vertebra, and one corresponding to the lower endplate 
of the immediately above vertebra — there should be a peak in the minimum curvature projection 
curve, corresponding to the intervertebral disc. Therefore, endplate positions are sought which 
should satisfy these conditions, but at the same time meet all constraints about vertebra shape. 
The constraints about vertebra shape used in this invention are: 

a) The heights of vertebrae are within a given range; 

b) The height ratio of neighboring vertebrae is within a certain range; 

c) The ratio of the thickness of the inter- vertebral disc to the vertebra 
height is within a certain range; 

d) The angle between neighboring endplates should not exceed a 
certain value. 

e) The vertebra height variation should satisfy a global model. 

Items a), b) and c) have been already disclosed in 00P7820US, while items d) and e) are 
new to this invention. These constraints correspond to box 120 of Figure 1. 

Referring again to Figure 2, ranges of vertebra heights can be computed as being 
proportional to the vertebra widths, which can in turn be measured from the distance between the 
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left and right spine boundaries. Suppose the vertebra positions on the spine axis are denoted by 
e x e 2 , e N , where N is the number of endplates. Suppose the order of the endplate positions 
are from the lower part to the upper part of the spine body, and e x is the lower endplate of the 
starting vertebra in the detection, and e N is the upper endplate of the ending vertebra. Then the 
height of each vertebra is h x =e 2 -e,, h 2 =e 4 -e 3 , h v = e N -e N _ x with the number of 
vertebrae being V -N 12. The thickness of the intervertebral disc is then 
h =^3 ~ e 2> f i = e 4 _e 3> •••> h-x =e N-i ~ e N-2- Suppose the orientation of endplate e v is 
denoted by 9 (e v ). Then the conditions a) through e) can be expressed, respectively, as 

Km <K < h ™ (!) 

(2) 

r th,r™ <t v l K <r th,™ ( 3 ) 

|^k)-^i)|<^ ( 4 ) 
g(K)=o (5) 



where , h m , r Mmin , r Mmax r th ^,r th ^ , and 9 m are predefined values. Because there exits 
no standard form for g () in physiology for describing the relationships among vertebra heights, 
we found the following generic form of g () gives a good approximation of this relationship 



h v =ka v +b 



(6) 



where the parameters k,a,b are unknown constants for a specific spine, but varies across 
different spines. Because of the order of endplates as being from bottom to top, the parameter 
a will always be less than or equal to 1, meaning that running from the lower to the top part of a 
spine, the vertebra height should never be increasing. If a = 1 , Equation (6) models vertebrae of 
5 constant heights. 

Based on the above analysis, the task of endplate detection can be formulated as: to find 
e x e 2 e N such that a score measuring the degree of support for these positions to be on 

endplates be maximized while satisfying the constraints of Equations (1) through (5). We 
choose the score as: 

U 

m 

w 

S E = Z M + £ / max {C J>min (p) } (7) 

m 

O where V has the same meaning as before, i.e., the number of vertebrae. The first term in 

Li. 

Equation (7) is the sum of the values of the maximum curvature projections at the endplate 
p*5 positions, while the second term is the sum of the values of the minimum curvature projections at 
the intervertebral disc positions. Because the exact positions of the intervertebral disc are not 
important to the endplate detection, the second term of Equation (7) is formulated to say that the 
intervertebral disc should be constrained to lie between an upper endplate and the lower endplate 
of the immediately above vertebra and to have a peak in the minimum curvature projection. 
20 The maximization of Equation (7) subject to the constraints Equations (1)~(5) is a 

nontrivial task. First, because the problem is highly nonlinear, there is no closed-form solution 
for it. Secondly because the dimension of the problem is N + 1 , an exhaustive search is 
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impractical. Note that the number of endplates N is also an unknown (that's why the dimension 
of the problem is N + 1 instead of N ). Finally the global constraint of Equation (6) with 
unknown model parameters complicated the problem further. 

We utilize a generalized dynamic programming approach for the optimization. In the 
standard dynamic programming approach, constraints like Equation (6) cannot be enforced. 
Because the model parameters of Equation (6) are unknown, they have to be estimated from the 
heights of all vertebrae. This means that all endplate positions are involved in Equation (6), 
making the constraint an N -dimensional one. In this invention, an approach is proposed to 
introduce backward tracing during forward propagation, so that global constraints like Equation 
(5) can be enforced in a progressive manner. 

First Equation (7) needs to be converted into a form suitable to solve by dynamic 
programming: 

+ Q max (e 2 )+C Sjmax (e,) + max {C sMn (p)} 

max 

(e 5 )+ max {C min (»}+... r 8 ) 
+ C s , I[m (e N _ 2 )+C s>aM (e N _ l )+ max {C sMn (p)} 

V 

= C ,,max (*1 ) + S f C ^max fev ) + Qmax fev+l ) + ™aX {C,^ ( P )}] + C s ^ x (e N ) 
v=l pe(e 2v e lv+l ) 

where the score terms are grouped in V + 2 groups. The maximization of Equation (8) can be 
viewed as choosing appropriate variables from V + 2 layers, e x from the first layer, 
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( e 2v> e 2v+i )from the following V layers, v = l 9 ...V , and e N from the last layer, such that the total 
score is maximized. We labeled the layers by Z v , v = 0,1,.. F + lin the sequel. 

Instead of directly using endplate positions as the state variables in each layer (except for 
layers Z 0 and I v+1 ), we employ an equivalent representation (e 2v ,t v ), where t v is the thickness of 
the intervertebral disc between vertebrae vand v + 1 . The advantages of using (e 2v ,t v )to 
parameterize the state space is that the range of t v is much smaller than that of e 2v+1 , reducing 
significantly the memory requirement in the maximization procedure. According to Equation 
(3), the range of t v is only a fraction of the vertebra height. From (e 2y 9 t v ), it is easy to get the 
other endplate position as 

e 2 v + i =e 2v +t v (9) 
With the new parameterization, the local score for layer L v in Equation (8) can be written 

as: 



EoM^C^M (10) 

pe(e 2v ,e 2v +t v ) 



Ev +l {e r )=C Stiaa (e v ) (12) 



Because at most two variables appear in each local score, the dimension of the dynamic 
programming is two if there were no constraint like Equation (6) involved. As was mentioned 
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before, the constraint of Equation (6) involves all N endplate positions, making the dimension of 
dynamic programming the same as that of an exhaustive search. 

The score maximization of Equation (8) by dynamic programming without the 
constraints of Equations (1)~(5) can be expressed as follows: 

S x {e 2 ,t x ) = max[£ 0 (e x ) + E x (e 2 , t x )] (13) 
S v (e 2v ,t v )= maxfSU (e 2(v _ 1} , t v _ x )+ E v (e 2v , t v )] ;v-7, . V 

e 2 ,t v 

(14) 



^ ( e N ) = max[S r (e 2V ,t v )+E v (e N )] (15) 



e 2V ,t v 



max E = max S F {e N ) (16) 



Referring to Figure 5, this procedure will be called forward propagation in the sequel. 
The solution for endplate positions is obtained by back-tracing from the maximum of the last 
layer. Figure 5 graphically depicts how this recursive score computation procedure can be 
viewed as an optimal path-finding problem, optimal in the sense of maximum score. 

Referring to Figure 6, in the presence of the constraints of Equations (1)~(5), we modify 
the forward propagation as follows. First, local constraints of Equations (1)~(4) are used to 
define the parameter ranges over which the maximizations are made. Here it is assumed that the 
ranges of the first and the last endplates e x and e N are known a priori at control box 610. 

Denote these ranges by [e l niin e l imx \ and [e N mhl e N ^ x J, respectively. (Note that the total number 

of endplates, N , is unknown.). From the height constraint (1), the range of e 2 can be obtained as 
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kmin, ^maxj=kmm + ^rmn ^l,max + A max J- For each C 2 € J, and each 

e i 6 k ; min ><W J > the corresponding vertebra height is ^ = e 2 - e x . Based on Equation (3), the 
range [t l mm 9 t hm ] for t x can be computed at control box 620. At control box 630, the propagation 
from the range [e vjn ^ e vmax J x [t vmin t VfWaL j at layer L v to the layer L v+l is obtained similarly. 

At control box 650, for the global constraint of Equation (6), local backward tracing is 
introduced during the forward propagation. For the maximization at layer L v , back-tracing is 

done for each element e [e 2{v ^ ,e 2(v _ lW Jx [t v ^ ,* v _ lmax j at layer . This identifies a 
locally optimal path for each p v _ x . From the local path of p y _ x , a profile of vertebra heights 
Hpv-i = {hpv-i (0> z = 1A.., v- 1} can be obtained. Then, at box 660, for each element 
/? v = (e 2v 9 t v ) at layer L v , we check whether the augmented height profile = {H^ , (v)} , 
with A^, (v) = e 2v - (e^j + ^ ), satisfies the global model of Equation (6). To do this, an 
iterative, alternating estimation procedure is adopted to estimate the parameters k, a, b first, 
based on the given profile. Initially, the parameter a is set as a (o) = 1 , where the super script 
represents the iteration number. Then, Equation (6) becomes linear in k and b , and can be 
solved for k {l \ b {l) . Inserting the obtained k {l \ b {l) into Equation (6), Equation (6) can be solved 

for a to get a a (l) . With the obtained new values £ (2) and 6 (2) are solved for. This 
procedure repeats until convergence is reached. With the obtained fitting parameter, the residual 
error for each height is checked at box 670. A height profile h{i\ i = l v .. ? K is said to satisfy the 
global model if Vz 

\h(i)-ka l -b\/h(i)<S h (17) 
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where S h is a predefined value. If the height profile satisfies the global model, element p v \s 
kept. Otherwise, it is discarded. All elements that remain at layer L v satisfy a global model and 

are then used at box 680 in the forward score propagation at the next layer. In this way, the 
global constraint is ensured in a progressive way. 

When the last endplate is scored, control branches at box 640 to box 645 where the 
maximum score is found for the last layer. At box 655, we backtrack from the position of this 
maximum score to find the endplate. No further operation is performed on the height profile, 
because it has already satisfied the model. 

Figure 7 shows an example of endplate detection using the disclosed technique in this 
invention. 

The detected endplates can be used to compute the wedge angle and tilt angel of any 
vertebra as well as the Cobb angle. All these angles are important measurements in the 
deformity analysis of spines, upon which physicians rely to make decision about diagnosis or 
surgical planning. 

The methods of the invention may be implemented as a program of instructions, readable 
and executable by machine such as a computer, and tangibly embodied and stored upon a 
machine-readable medium such as a computer memory device. An embodiment may be in the 
form of a computer controlled module into which is input a spinal image such as Figure 7(a) and 
returns the image superimposed with the endplate positions, such as Figure 7(b). 

It is to be understood that all physical quantities disclosed herein, unless explicitly 
indicated otherwise, are not to be construed as exactly equal to the quantity disclosed, but rather 
as about equal to the quantity disclosed. Further, the mere absence of a qualifier such as "about" 
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or the like, is not to be construed as an explicit indication that any such disclosed physical 
quantity is an exact quantity, irrespective of whether such qualifiers are used with respect to any 
other physical quantities disclosed herein. 

While preferred embodiments have been shown and described, various modifications and 
substitutions may be made thereto without departing from the spirit and scope of the invention. 
Accordingly, it is to be understood that the present invention has been described by way of 
illustration only, and such illustrations and embodiments as have been disclosed herein are not to 
be construed as limiting to the claims. 
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